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ABSTRACT 

Tempo2 is a new software package for the analysis of pulsar pulse times of arrival. In 
this paper we describe in detail the timing model used by tempo2, and discuss lim- 
itations on the attainable precision. In addition to the intrinsic slow-down behaviour 
of the pulsar, tempo2 accounts for the effects of a binary orbital motion, the secu- 
lar motion of the pulsar or binary system, interstellar, Solar system and ionospheric 
dispersion, observatory motion (including Earth rotation, precession, nutation, polar 
motion and orbital motion), tropospheric propagation delay, and gravitational time 
dilation due to binary companions and Solar system bodies. We believe the timing 
model is accurate in its description of predictable systematic timing effects to better 
than one nanosecond, except in the case of relativistic binary systems where further 
theoretical development is needed. The largest remaining sources of potential error 
are measurement error, interstellar scattering, Solar system ephemeris errors, atomic 
clock instability and gravitational waves. 
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1 INTRODUCTION 

The core of any pulsar timing software package is the mathematical model for the times of arrival (TOAs) of pulses at 
the observatory. As the pulsar population grows more diverse and observing systems increase in sensitivity, the accuracy 
demanded of timing models grows accordingly. For many years the standard choice for pulsar timing was the tempo 1 software 
pac kage. Accurate at the 100-ns level, the timing model of tempo was adequate for numerous important disco veries (e.g. 
ITavlor et ailll979l: iTavlor fc Weisber j Il989t IWolszczan fc Fraill Il993 : iKaspi et al1ll994l : Ivan Straten et al.ll200lTh However, 
motivated primarily by the possibility of directly detecting gravitational radiation through its effect on pulse arrival times, 
current observing campaigns are setting highly ambitious goals in terms of data precision and volume. These projects aim to 
obta in of the order of ~ 10 4 independent TOA measurements with a typical root-mean-square (RMS) uncertainty of 100 ns 
(e.g. Ijenet et al.ll2005l) . dictating a maximum allowable systematic error of 1 ns. 

Motivated by the need for improved accuracy as well as greater flexibility in fitting and analysis procedures, we have 
developed a new p ulsar timing package called tempo2. The architecture, capabilities and usage of tempo2 are described by 
Irlobbs et alJ J2Q06) (hereafter Paper I). The purpose of the present work is to describe the timing model in detail, and provide 
estimates of the accuracy of each component of the model. Section [5] begins with an overview of the model and definitions of 
the various relativistic coordinate frames and proper times (Section 12.21 . This is followed by an expansion of the complete 
geometric propagation delay accurate to better than 1 ns (Section 12.3^ . The remainder of Section |5] discusses all oWther 
delay terms in the timing formula, and provides a detailed description of how each contribution is evaluated according to 
the parameters of the model and other information. Section El gives a detailed account of the accuracy with which each term 
is computed, and discusses possible avenues for further improvement of the timing model. Finally, for the convenience of 
comparison with the commonly used tempo software package, Section [1] details the differences between tempo and tempo2. 
Tabular summaries of timing model parameters and variables are provided in the Appendix. 



* Email: Russell.Edwards@csiro.au 

1 http:/ /www. atnf.csiro.au/research/pulsar/tempo/ 
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2 TIMING MODEL 

2.1 Overview 

Pulsar timing consists of the measurement of the pulse times of arrival (TOAs) at the observatory and the fitting of these 
TOAs to a model. The timing model relates the measured TOA to the time of emission at the pulsar, from which a "pulse 
phase of emission" is computed via a model of the intrinsic variations in the pulse period. If the timing model is correct, 
in the absence of measurement noise all calculated pulse phases will be integer (when expressed in units of cycles or turns). 
In practise, the parameters of the model are not known a priori, and measurement noise is present, so that the deviations 
of these pulse phases from the nearest integer (the "residuals") are finite. Once the timing model parameters are refined to 
optimise some statistic of the residuals - in tempo2 this is the root-mean-square or \ 2 ( see Paper I) - and assuming there 
are no shortcomings in the structure of the timing model, the residual deviations are due only to measurement noise. 

The focus of this paper is the timing model, which relates a measured TOA to a time of emission, and from there computes 
a corresponding pulse phase of emission. The latter step is described in Section [2.81 The former part comprises the bulk of 
the complexity of the model, and it too can be helpfully broken into two main parts. Firstly, because that part of the model 
relates a proper time of arrival at the observatory to a proper time of emission at the pulsar, a general relativistic frame 
transformation is necessary. This is performed in a sequence of manageable steps, outlined in Section f2. 21 Secondly, the light 
travel time must be accounted for. This is derived in Section 12.31 under the assumption of Euclidean space-time and zero 
refraction or dispersion. The remaining parts of Section [5] derive in full mathematical detail the various parts of the frame 
transformation and light travel time, including dispersion, refraction and general relativistic effects along the propagation 
path. 

2.2 Coordinate frames and time variables 

The relativistic frame transformation between observatory proper time and pulsar proper time is computed in several steps. 
This section outlines the steps, gives the names of the various reference frames and introduces the notation use to refer to 
time variables in each frame. The time variable for each frame may appear subscripted with the letter "a" or "e", where "a" 
denotes a time of arrival and "e" a time of emission, at the origin of the frame denoted by the superscript. A time variable may 
also appear without a subscript, meaning that it is an an arbitrary coordinate time variable. As such, time variables without 
subscripts are related by coordinate transformations, whereas subscripted time variables additionally include propagation 
delays. Spatial 3- vectors are defined in terms of a right-handed Cartesian system. 

The pulsar timing procedure begins with a measured arrival time, which after correcting for imperfections in the obser- 
vatory clock yields a time of arrival t° bs . The clock correction process in tempo2 operates through interpolation of tables of 
measured offsets between pairs of clocks and their variation with time (see Paper I), and includes offsets relating to the fact 
that observatory clocks usually approximate Coordinated Universal Time (UTC 2 ) rather than Terrestrial Time (TT). The 
corrected arrival time is measured in the TT time scale, which is defined to differ from coordinate time at the geocentre by 
a constant rate. It should be noted then that times referred to TT are only approximately proper times of the observatory. 
This transformation to Geocentric Coordinate Time scale (TCG) is done in the same step as transformation to Barycentric 
Coordinate Time (TCB), so no variable is assigned to this time scale. The spatial part of the terrestrial frame is the Inter- 
national Terrestrial Reference System (ITRS), while the geocentric counterpart is the Geocentric Celestial Reference System 
(GCRS). 

Coordinate times t SSB in the frame of the Solar system barycentre (SSB) are obtained by application of the gravitational 
redshift and special relativistic time dilation integral ("Einstein delay") described in Section fe.^.ll Coordinate time (TCB) in 
this frame represents the proper time that would be measured by an observer located at the Solar system barycentre were the 
gravitational field of the sun and planets not present. The spatial part of the barycentric frame is the Barycentric Celestial 
Reference System (BCRS), while directions in this frame define the International Celestial Reference System (ICRS). The 
origin of the BCRS is the SSB, hence if SB refers to the arrival time of the pulse at the barycentre. 

For further information on the definitions of TT, TCG, TCB, the ITRS, GCRS, BCRS and ICRS, see the Resolutions o f 
the XXIIIrd and XXIVth General Assembly of the International Astronomical Union (IAU; lAndersorJll99flt lRickman|l200lh 
and the 2003 Conventions of the International Earth Rotation and Reference Systems Service ttER.Si lMcCarthv fc Petitl200^ . 

Tracing the propagation backward toward the source, the next reference frame is that which is co-moving with the 
pulsar, or if it is a member of a binary system, the binary barycentre (BB). Coordinates in this frame are defined such that 
the coordinate time t BB would be equal to the proper time of an observer at the origin of the frame, were the gravitational 
field of the pulsar and its companion not present. The transformation of the arrival time to this frame is discussed in Section 

2 Due to copyediting errors, in Paper I the acronyms UTC, TCG, TCB, and TDB were inadvertanty referred to as UCT, CGT, BCT 
and BDT. 



Tempo2: Timing model and precision 3 



12.6.31 The orientation of this frame is chosen such that the Lorentz transformation relating it to the SSB frame is rotation-free, 
so that to zeroth order the corresponding basis vectors are parallel. The spatial origin of this frame is the BB. The time origin 
is chosen such that t BB has zero offset from t SSB at t BB = t SSB = t pos , where the latter is an epoch of position set by the user. 

The time scale for measurements is the proper time t psr measured at the pulsar, or more precisely at its centre, were 
its gravitational field not present. For isolated pulsars t psr = t BB . For binary pulsars, t pBT — t BB is the binary Einstein delay 
(Section 12.7.40 . As with the SSB-BB case, the BB-pulsar transformation is rotation-free. The spatial origin of this frame is 
the centre of gravitational field of the pulsar, so that t psr is the equivalent time of emission from the centre of the star. This 
differs from the actual time of emission by (to first order) the light travel time to the true point of emission, which is assumed 
to be stable on the time scale of interest. 



2.3 Geometric propagation delay 

Although it is convenient to break the propagation delay into separate contributions and treat them in the order in which 
they occur along the photon path, some terms in the model depend jointly upon the motion of the observatory and the motion 
of the pulsar (relative to the SSB). For this reason, in this section we derive the complete geometric propagation delay, that 
is the time of propagation neglecting dispersion, refraction and relativistic delays. The resultant terms are then assigned to 
different parts of the propagation path and treated in the relevant sections below. 

The geometric propagation delay is simply the Euclidean distance from the observatory to the pulsar, divided by the 
vacuum speed of light, c. The displacement vector from the observatory to the pulsar is the sum of the bary centric position 
of the observatory (r), the barycentric position of the pulsar (or BB) (Ro) at a given epoch (t pos ), the position of the pulsar 
with respect to the BB (b, zero for isolated pulsars), and the displacement (k) of the binary barycentre (or the pulsar itself, 
if isolated) in the time elapsed since the epoch t pos , owing to the initial velocity of the system (relative to the SSB) and 
acceleration in the Galactic or cluster gravitational field: 

R = R + b + k r. (1) 
Squaring equation Q and expanding the square root to third order, 

\R\ = [\Rof + \b\ 2 + \k\ 2 + \r\ 2 + 2(R -b + R -k- R -r + b- k -b-r -k -r)] 1/2 (2) 
A 2 A 3 

\ R °\+ A -^-l + ^-H + ---> where ( 3 ) 
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-i-r \\b\ 2 + |fe| 2 + |r| 2 + 2 (Ro ■ b + R ■ k - R ■ r + b ■ k - b ■ r - k ■ r)] . (4) 
Ro\ 



Neglecting terms of order \Ro\ 3 and denoting the initially radial and transverse components of each vector using subscripts, 
i.e. aii = a • Ro, a± — a — a^Ro/\Ro\ and hence, oi ■ b± — (a X Ro) ■ (b X Ro), the following relation is found: 

\R\ = |.Ro] +fc|| -r|| +6|| 

+ra (¥ — + ^ ^ ¥ - ■ * + ¥) (• - & + & - &) + °<'"-r 3 >' <«> 

The first four terms in equation JSJ are the initial SSB-BB distance (Section 12.6.21 . the secular displacement in the initially 
radial direction (Section 12. (i. 21 . the projection of the observatory-SSB vector on the initial line of sight (Section 12.5.31 . and 
the projection of the binary motion on the initial line of sight (Section 12.7.11 . The terms in the first pair of parentheses 
correspond to the Shklovskii effect 12.6.21 . annual proper motion 12.5.31 . secular changes in the apparent orbital viewing 
geometry 12.7.11 . annual parallax 12.5.31 . annual-orbital parallax 12.7.11 and orbital parallax 12.7.11 . The second, third and 
fourth terms in the second pair of parentheses represent variations in the effects just discussed due to radial motion of the 
pulsar or binary system, the motion of the Earth around the SSB, and the motion of the pulsar about the BB. Although 
these terms give rise to an additional eighteen contributions to the timing formula, to reach the stated accuracy goal of 1 ns, 
for a 20-yr observing campaign on any of the presently known millisecond pulsars only three of them need be considered: the 
secular change (Section 12.6.21 and annual modulation (Section 12.5.31 of the Shklovskii effect and the secular change of the 
annual proper motion (Section l2.5.31 . 

The displacement due to secular motion, k may be broken into first and second derivatives: 

k = n\Ro\ (t BB - tpo S ) + ^ (t BB - tpos) 2 , (6) 

where fi is the velocity divided by the distance, or in essence a three-dimensional proper motion, and a is an acceleration 
vector accounting fo r Galactic differential rotation and gravitational acceleration. The acceleration is of order ~ 10~ n m s~ 2 
jBell fc Bailed Il996h . contributing only ~ 500 km to h at either end of a 20- year observation campaign centred on £ p0 s- 
Equation JS) can be expanded in terms of this parameterisation of k, giving rise to an additional 27 terms. We have estimated 
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the magnitude of these terms for all pulsars in the ATNF catalogue 3 (Manches ter et al J 12005*1. including an allowance for 
greater acceleration for pulsars in globular clusters (a ~ 10 m s" " 2 ; lFreire et alJbOQll) . Only two of the extra terms involving 
a exceed 1 ns for a 20-year observing campaign, the first contributing to fcii and the second contributing to the Shklovskii 
term of equation ©. For all other terms of equation 10 involving k the acceleration is neglected. 



2.4 Top-level timing formula 

Having described the coordinate frames used and derived the basic Euclidean propagation length, it remains to specify the 
details of the calculation of the various delay terms associated with frame transformations, the vacuum propagation delay, 
and corrections due to non-unity refractive indices and space-time curvature along the propagation path. For convenience the 
timing formula is divided into three main contributions, which will be further broken down in the sections that follow: 

tT =tf s -A -A IS -A B . (7) 

A© f Section 12. 5H includes the coordinate transformation to the SSB frame, vacuum propagation delays associated with 
the Earth's orbital motion and the spin, precession and nutation of the Earth (on which the observatory is located) and 
excess delays owing to the passage of the signal through the Earth's atmosphere and the Solar system. This term relates the 
measured time of arrival to the TCB arrival time at the SSB: r| SB = t° hs — A Q . 

Ais 12.61 1 includes the transformation to the binary barycentre frame, vacuum propagation delays due to the secular 
motion of the system, and excess propagation delays due to passage of the signal through the interstellar medium. This term 
relates the TCB time of arrival at the SSB to the BB coordinate time of arrival at the BB: t BB — t^ SB — Ais. 

Ab includes the transformation to the pulsar frame, vacuum delays due to the binary orbital motion, and excess delays 
due to passage of the signal through the gravitational field of the companion. This term relates the BB coordinate time of 
arrival at the BB to the pulsar proper time of emission: tg BT = t a B — Ab- 

Note that since some of the geometric delay terms involve two of the three components (Solar system, secular or binary) 
of the relative motion of the observatory and the pulsar, the choice of which of these terms are assigned to each of the three 
top-level terms is somewhat arbitrary. See Section f2.3l for a detailed listing of where each term is discussed. 

Once the time of emission is obtained, it remains to compute the phase of the pulse train at that time. This relationship, 
4>(tc ST ), is discussed in Section 12.81 For the correct set of model parameters and in the absence of measurement noise, the 
phase for every computed time of emission is an integer number of cycles (turns). The fractional part of <p(t) is the timing 
residual, which is more often expressed in time units through division by the pulse frequency. 



2.5 Forming bary centric arrival times 

The term denoted A© in equation Q includes all effects pertinent to the transformation of the observed TOA (t^ 35 ) to an 
equivalent time of arrival for the same pulse wave-front at the Solar system barycentre. This in turn can be broken into several 
steps. These are: atmospheric delays, vacuum retardation due to observatory motion ("Roemer delay" and parallax), excess 
propagation delay due to dispersion, the effects of relativistic frame transformation ("Einstein delay"), and the excess delay 
experienced by rays as they pass through the gravitational potential of Solar system bodies ("Shapiro delay"). These terms 
are denoted as follows: 

A© = A A + A R0 + A p + A D © + A E © + As© (8) 



2.5.1 Einstein delay 

The space-time coordinates of pulse arrival events are specified in the coordinate frame of the observatory, in which pulse 
timing anomalies may manifest due to spin and orbital acceleration and variations in gravitational potential. To avoid these 
effects, pulse arrival events are transformed to the quasi-inertial frame of the Solar system barycentre. The two frames are 
related by a relativistic four-dimensional space-time transformation. The spatial part of the event is the displacement of the 
observatory from the geocentre at the instant of reception, which is altered by a negligible amount due to special and general 
relativistic length contraction. The eff ects of relativistic time dila tion, on the other hand, cannot be neglected. te mpo2 uses 
the n umerical time dilation results of llrwin fc Fukushimal dl99Sl) . who used the DE405 Solar system ephemeris istandishl 



1998) to compute the time dilation integral: 



A E ©_© = 1 / f U e + '"-f + AL<T N > + A4 A > ) At, , 9) 



3 http:/ /www. atnf.csiro.au/research/pulsar/psrcat/ 
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where [7® is the gravitational potential at the geocentre due to all Solar system bodies except the Earth, and is the velocity 
of the geocentre relative to the Solar system barycentre. The first two terms describe, to order 1/c 2 , the gravitational redshift 
and special relativistic time dilation respectively. These contribute in roug hly 2:1 ratio (±[/ e ~ v% ~ l(T 8 c 2 ), for a mean 
drift of ~ 1.5 x 10" 8 , or ~ 0.5 s yr" 1 . The main source of time variation in the integrand is due to the orbital acceleration 
of the Earth, which via conservation of energy causes approximately equal changes (of the same sign) in the first two terms, 
amounting to ~ 3 x 1CF 10 in rate amplitude, or ~ 2 ms in int egrated delay am plitude. The remaining terms apply a correction 
for h igher-order relativistic terms (AL^ N) = 1.097 x 10~ 16 : lFukushimalfl995]) and asteroids (AL^ = 5 x 10~ 18 : lFukushimal 
Il995l) . made in mean rate only. 

The time dilation integral of equation © relates the coordinate times of clocks at the geocentre and the Solar system 
barycentre. For observers located on the surface of the Earth, corrections must be made for the differential time dilation and 
gravitational redshift between the observatory and the geocentre, yielding: 

where s is a vector from the geocentre to the observatory, r® is the velocity of the geocentre with respect to the barycentre, 
and Wo = 6.969290134 x 10 _10 c 2 approximates the gravitational plus spin potential of the Earth at the geoid. (See Section 
I2.5.3l for the calculation of s and rm .) N ote that, while the value of the potential varies geographically and temporally, the most 
recent definition of TT (Rick marj200lL the terrestrial time scale to which tempo2 refers pulse time of arrival measurements, 
takes the exact value of Wo/c 2 quoted above as definitive of the rate difference between TT and SI coordinate time at the 
geocentre. Given a common epoch of to = 43144.0003725 (Modified Julian Date), equations 15*1 and 1101 relate measured TT to 
the IAU recommended coordinate timescales TCG ( geocentric) and TCB (bary centric) . 

It should be noted that the time ephemeris of llrwin fc Fukushimal II19991) in fact takes its time argument, T ep h, in a 
reference frame that is a linear transformation of TCB, making the computation of the Einstein delay strictly an inverse 
problem. In practise, since the argument time scale has zero mean rate difference to TT, the latter can be substituted with 
negligible error in the result lllrwin fc Fukushimal 119 991 . In addition, the evaluation of equation 1101 requires consulting 
the Solar system ephemeris (Section 12.5.31 which also takes T cp h times as its argument. Since T ep h is computed from TCB 
(equation llOH . a circular dependency problem arises. This is solved to sufficient accuracy by two passes through an iterative 
refinement process. 



2.5.2 Atmospheric Delays 

The group velocity of radio waves in the atmosphere differs from the vacuum speed of light. Refractivity is induced both by 
the ionised fraction of the atmosphere (mainly in the ionosphere) and the neutral fraction (mainly in the troposphere). 

The tropospheric propagation delay be separated into the so-called "hydrostatic" and "wet" components. In Very Long 
Baseline Interferometry (VLBI) and Global Positioning System (GPS) applications, these components are typically modelled 
as the product of the delay induced at the zenith, and a so-called "mapping function", which specifies the ratio between the 
zenith delay and the delay in a given direction. For a planar atmosphere, the mapping function is simply given by esc O where 
O is the source elevation angle. More advanced mapping functions take into account the curvature of the atmosphere, assuming 
azimutha l symmetry. The effect of azimuthal asymmetry is typically less than a nanosecond even at very low elevation angles 
(e.g. lMacMillanlll995l) . and need not be considered for the present purpose. 

The hydrostatic component contributes approximately 90% of the total d elay, and may be co mputed a priori assuming a 
given mixture of gases in hydrostatic equilibrium. tempo2 uses the formula of iDavis et alJ Jl985l) . re-written as a time delay 

A (43.921 kPa) / -1 -1 \ 

hZ ~ c [1 - 0.00266 cos 2y- ( Mg^jm )] ' 1 > 

where Ah z is the hydrostatic zenith delay, P is the surface atmospheric pressure, <p is the geodetic latitude of the site, and H is 
its height above the geoid. In combination with the mapping function, this leads to a timing term of amplitude ~ 7.7 ns ■ esc O, 
mostly on a diurnal timescale. If atmospheric pressure data are unavailable, tempo2 u ses a canonical value of one standard 
atmosphere (101.325 kPa). tempo2 uses the Niell Mapping Function (NMF: lNielllll996f) . which is of comparable accuracy to 
other published mapping functions but does not require meteorological data: it depends only on the source elevation. 

The wet component of the tropospheric propagation delay is highly variable and cannot be predicted accurately. Fortu- 
nately, it is small. The zenith wet delay (ZW D; A W z) may be m easured using appropriate analyses of radiosonde, water vapour 
radiometer, GPS or VLBI observations (e.g. iNiell et al.lf200lh . of which only GPS is likely to be obtainable on a routine basis 
at most radio observatories. If such measurements are available, tempo2 can make use of them in conjunction with the NMF 
to account for the effects of tropospheric water vapour on pulse times of arrival. The ZWD could conceivably be included as 
a free parameter in the pulsar timing model, but because the effect is small and varies from day to day, this would not be 
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possible with the sensitivity of presently attainable observing systems. For this reason, if no tabulated ZWD information is 
available, the effect is neglected. 

Neglecting dispersion (below), the total atmospheric delay is written: 

A A = m h (e)A hz (P)+m m (e)A wz (P), (12) 

where mj,(0) and m w (Q) are the Niell mapping functions for the hydrostatic and wet components. Values for P and A wz are 
to be provided as input data to tempo2, or default values of P = 101.325 kPa and A wz = are used. 

The passage of the signal through the ionosphere induces a dispersive delay that varies strongly with the Solar activity 
cycle and also on on seasonal, diurnal and shorter timescales. The integrated column density of electrons ("total electron 
content"; TEC) in the ionosphere typically lies in the range 5-100 TECU (1 TECU = 10 16 m~ 2 ~ 3.2 x 10" 7 cm -3 pc), 
corresponding to variable propagation delay s of the orde r of 7—130 ns -(//l GHz) -2 , where / is the observing frequency. 
Models of th e ionosphere ar e available (e.g. ISchaeJll99St iBilitzalbOQlh . however the predicted TEC is accurate to only ~ 
3-20 TECU (|Mamoru et aEH jjjj). In any case, ionospheric dispersion is inseparable from and smaller than uncertainties 
in interplanetary (Section I2.5.4|l and interstellar (Section 12.6. 1^ dispersion, necessitating the inclusion of a fittable, fully 
time- variable dispersion parameter ('Section l2.6.1^ . 



2.5.3 Roemer delay and parallax 

The Roemer delay is the simple vacuum delay between the arrival of the pulse at the observatory and the Solar system 
barycentre, not including effects related to the binary motion of, or finite distance to, the pulsar: 

A R0 = -^K (13) 

where Rbb = -Rbb/|-Rbb| is a unit vector in the direction of the binary barycentre at the time of observation. tempo2 
performs this calculation in the BCRS. 

The vector Abb is constructed from the spherical coordinate angles of the ICRS source direction, right ascension (a) and 
declination (8), at time t pos , and the Cartesian components of the proper motion of the source in the plane of the sky (fi a 
and jj,s) and along the line of sight (/xii), all of which are fittable parameters in tempo2: 

Rbb = -Ro + M-L (*a B - *pos) - Q|M-l| 2 -Ro + (*a B - *pos) 2 ,where (14) 

cos a cos 8 

Ro = | sin a cos 8 | (15) 
sin 8 

(J-± = Mq« + M<5<5, (16) 
and (17) 

- cos a sin 8 

- sin ol sin 8 | . (18) 
cos 8 

Alternatively, ecliptic coordinates (A, (3, may be used, in which case: 

cos A cos (3 

Ro = E | sin A cos (3 | , and (19) 

sin /3 

M-L = MaA + fjLp/3, where (20) 
(X \ 

E = coseo — sineo , (21) 






A = E [ cosA j , and (22) 
- cos A sin (3 

= E | -sinAsin/3 | . (23) 

cos (3 
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Here eo = 84381.40578 arcsec faarada fc Fukushimal l2004"l is the mean obliquity of the ecliptic at J2000.0. 

The four terms that result from the expansion of equation 1141 derive from corresponding terms in the expansion of 
equation JSJ, such that 

- k± r± nijfcj-l 2 k«k±-r± ,„.. 

Here we have assumed that k = /x|J?o| (t BB — ipos), safely neglecting any acceleration or higher order secular motion of the 
pulsar or BB ( Section 12. 3^ . In principle, t BB cannot be computed without first knowing Arq; this is resolved by initially 
setting t BB to the TCB time of arrival at the observatory, and iteratively refining Aq and Ais until convergence. 

Equation 1241 consists of four terms corresponding to the expression in equation 1141 of the instantaneous source direction 
as the sum of four contributions. These are the initial direction (yielding the first term of equation 1241 . the proper motion 
(second term), a correction to maintain unit length after addition of the proper motion (third term) and the acceleration or 
deceleration of the proper motion owing to the fact that over time, the proper motion moves the source vector, causing part 
of the initially radial velocity to acquire a transverse component (fourth term). The third term can also be interpreted as the 
propcr-motion-induced second-order reduction in the projection of the source direction on the radial component of r that 
accompanies the first-order increase in its projection upon the transverse component of r. Another alternative interpretation 
of the third term is that it corrects the Shklovskii term (Section 12.6.21 for the distance variations induced by the Earth's 
orbital motion. The fourth term, first introduced by N. Wex (1997, unpublished contribution to tempo), provides one of 
two possible ways of accessing the elusive radial velocity, the integrated first-order contribution of which (ku) is otherwise 
inseparable from a Doppler shift of the pulse and orbital frequencies and their derivatives. The second method is described in 
Section 

The first term of equation 1131 describes the difference in time of arrival of a pulse at the SSB compared to that at 
the point obtained by projecting the SSB-observatory vector on the SSB-BB vector. The arrival time at that point differs 
again from the arrival time at the observatory due to the curvature of spherical "wavefronts" connecting photons emitted 
simultaneously from the pulsar. This term is referred to as the "parallax" term, although it differs somewhat from stellar 
parallax: positional astrometry measures the three-dimensional orientation of the wavefront normal, whereas pulsar timing 
measures the position of the wavefront in its intersection with the ecliptic. To first order the curvature is proportional to the 
square of the lateral displacement of the observatory from the SSB-BB vector, yielding an excess delay corresponding to one 
of the terms of equation @: 

Here the "parallax distance", d p is substituted for \Ro\ to allow this effect to be separated from others involving the source 
distance. The conventional "parallax" n = 1 AU/d p , being the angle subtended from the source by line of 1 AU in length at 
the SSB, is a fittable parameter in tempo2. Its value is specified in units of milliarcseconds, corresponding numerically to the 
reciprocal of the distance in kiloparsecs. 

In computing the Roemer and parallax delays, the vector r is constructed in two steps, from the SSB to the geocentre and 
from the geocentre to the observatory, i.e. r — r® + s. The first part (rm) is accomplished with the aid of a numerical Solar 
system ephemeris. The current default ephemeris is DE405 l|standishlll998h . which is align ed within a kn own uncertainty to 
the ICRS, but uses a time and length scale slightly d ifferent than the SI second and metre JStandishlll998i) . tempo2 uses the 
equations and constants of llrwin fc Fukushimal Jl999h to appropriately transform the TCB site arrival time fSection l2.5.11 for 
input to the ephemeris (using an offset and scale) , and appropriately scales the output vectors to SI units. It should be noted 
that the use of alternative ephemerides that require a different transformation procedure will introduce errors (see Section 

Egj) . 

The second part of the observatory position vector, s, is obtained by transformation of the terrestrial observatory co- 
ordinates into the GCRS (which to sufficient accuracy is equivalent to the BCRS). This is accomplished using an ITRS 
site position vector (arms) supplied to tempo2, transformed to the GCRS using algorithms consistent with the IAU 2000 
Resolutions. The transformation proceeds as follows: 

a = Q(t° bs )R(e s )W(t° b > ITR s, (26) 

where the matrices W, R and Q account for polar motion, Earth rotation, and the motion of the Earth spin axis in 
the ICRS. The latter term in turn consists of a product of a frame bias and the IAU2000B precession- nutation matrix 
jM cCarthv fc Liizumfcoosl) . Details of the construction of these matrices are available in the IERS Conventions ^McCarthy fc Petitl 
2004) and in the documentation of the IAU Standards of Fundamental Astronomy (SOFA) software library 4 used by tempo2 
to compute them. Polar motion and Earth rotation both exhibit unpredictable variations that must be corrected post facto 



4 http://www.iau-sofa.rl.ac.uk/ 
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using measured parameters. tempo2 uses the "C04" series of Earth orientation parameters, provided by the IERS, as input 
to the SOFA routines. 

For the purposes of transforming the observing frequency to the barycentric frame fSection l2. 6. II . the time derivative of 
the Roemer delay is needed: 

dAfl _ (r 9 +s)-Rg 
dt ~ c ' 

where dots denote the time derivative and the effects of proper motion and parallax are negligible. The barycentric velocity 
of the geocentre, r®, is provided by the Solar system ephemeris. The barycentric velocity of the site is dominated by the 
velocity imparted by Earth rotation. This is calculated using 

sitrs = (wffiW _1 i) X SlTRS- (28) 

Here a unit North vector (z) is transformed using the inverse polar motion matrix and multiplied by the instantaneous Earth 
angular rotation rate to give the angular velocity of the Earth in the ITRS. The cross product of this with the site radius 
vector gives the tangential velocity, which is transformed to the GCRS in the same manner as s. Of the neglected terms in 
this approximation, the largest is the precessional velocity, at a negligible ~ 10 _13 c. 



2.5.4 Solar system dispersion 



Radio signals encounter significant dispersion in the interplanetary medium, due to the e lectron content of th e Solar wind. 
The electron distribution follows a roughly r~ 2 form consistent with spherical expansion Jlssautier et aljfl998l) . Integrating 
along the line of sight yields a dispersion measure (DM): 



no 



n (1 AU) 



1 AU 

r(s) 

2 P 



ds 



(29) 
(30) 



\r\ sin p 

where \r\ is the heliocentric radius of the observatory, p is the pulsar-Sun-observatory angle, and no is an overall scale parame- 
ter, which can be specified in tempo2 parameter files as NE1AU. (It should be noted that equation <■'-!( H is an approximation that 
ignores the effects on the dispersion relation due to the bulk plasma velocity and its variation along the line of sight. However, 
there are other much more significant sources of error as discussed below.) The defau lt value for no is 4 cm -3 , consistent with 
recent measurements Jlssautier et aljll99St llssautier et aljBoOll ISplaver et al.ll2005F) and significantly lower than the value of 
9.961 cm" 



iD0 



" 3 used by tempo. The computed dispersion measure is used to calculate the delay due to interplanetary dispersion: 
DM Q 



2.410 x 10- 16 cm- 3 pc 



(/ 



SSD\ 



(31) 



where / SSB is the observing frequency transformed to the barycentric frame (see Section \'2. 6. II . 

This correction is of limited accuracy. Recent studies of with the Ulysses spacecraft have revealed extraordinary complexity 
in the Solar wind. At Solar minimum, th e electron density (sc aled by r 2 ) was found to show RMS temporal variations of 10— 
50%, depending on heliocentric latitude l|lss i ajrtie i r - ^^uj^9 i 8]li while at Solar maximum the mean electron density increased 
but was subject to very strong modulation l|lssautier et all 2001). The default correction of equation l|31fl is therefore uncertain 
within at least factor of two, corresponding to a minimum of 130 ns of error at / = 1 GHz for a source located at an ecliptic 
pole, increasing to many microseconds for sources within a few degrees of the ecliptic. 

Errors in the predicted Solar system dispersion delay are of concern not only because they add random noise to the timing 
residuals and fitted parameters, but also because a porti on of the error term will exhibit annual periodicities that could corrupt 
the fitted pulsar position, proper motion and parallax (ISplaver ■ 



It is also p ossible that other pe riodicities in this 
due to a planetary companion to the pulsar ((Scherer et al.llT997L though see also 
,s to improve the model have focussed on adding fittable degrees of freedom (e.g. 
however the complex behaviour seen in the Ulysses observations indicate that the 
only way to adequately remove the effect is to measure it directly using multi- frequency observations (see Section \'2. 6, It . 



term could be misinterpreted as arising 


Wolszczan et al. 2000 


). Previous attemp 


Coenard et al. 1996: 


ISplaver et alJl2005h 



2.5.5 Shapiro delay 

The Shapiro delay accounts for the time delay caused by the passage of the pulse thro ugh curved spacetime. The total delay 
obtained by summing over all the bodies in the Solar system (backer fc Hellingslll986T) : 

A S0 = -2 Y, ^ ■ r i + \*i I) + A so2 (32) 

3 
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where R is a unit vector in the direction of the pulsar, rrij is the mass of body j, rj is a vector from body j to the telescope, 
and As02 is a second-order correction discussed below. The first term may be re-written in the following more convenient 
form: 

As ° = " ^r 1 ln |rj 1 (1 ~ cos ^ ) + As02 - (33) 
j 

where ipj is the pulsar-telescope-object angle for the j-th object. Tempo2 includes the first-order Shapiro delay for the Sun, 
Venus, Jupiter, Saturn, Uranus and Neptune. 

For observations of sources very close to Solar system bodies, higher-order effects may become impo rtant. The largest of 



these is the geometrical excess Dath length due to gravitational light bending (Richter & 


Matzncr 1983), which amounts to 


9.1 ns for a rav with a traiectorv grazing the Solar limb 5 . From the derivation oflHellingsl 


|l986|), assuming general relativity 


and re-parameterising in terms of ip for a very distant source, this term is given by 




. 4G 2 m 2 

A S02 ~~ si u ; • ; 

c°\r \ tan tp sm tp 


(34) 


4G 2 m 2 
c 5 \r\ip 2 


(35) 



For pulsar timing purposes, only the geometric delay due to the Sun need be considered. 



2.6 Propagation through interstellar space 

The propagation time of the signal from the pulsar (or binary system barycentre) to the SSB is conveniently broken into three 
contributions: the vacuum propagation delay, i.e. the path length divided by the vacuum speed of light, and the excess delay 
due to dispersion and other frequency-dependent delays: 

Ai S = A V p + Aisd + A PDD + A ES . (36) 

The fourth term accounts for the special relativistic time dilation between the SSB reference frame and the binary reference 
frame. 



2.6.1 Interstellar dispersion and other frequency-dependent delays 

The interstellar dispersion delay enters the timing model under the standard relation: 

(37) 



(/ SSB )2 > 

where D is the so-called dispersion constant, and / SSB is the frequency of the radiation at the Solar system barycentre. The 
barycentric frequency differs from the frequency at the observatory, due to a simple Doppler shift (~ 10~ 4 in magnitude) and 
a second-order relativistic correction (~ 1CP 8 ), plus the effects of gravitational redshift (~ 10~ 8 ). The first term is simply the 
time derivative of the Roemer delay (the "Roemer rate"), while the sum of the latter two are is derivative of the Einstein 
delay (the "Einstein rate"). 

/ SSB =(! + ^f + )/• (38) 

Although D is the directly measurable quantity of equation 1371 , it is customary to speak in terms of the inferred column 
density of electrons along the line of sight, i.e. the dispersion measure, DM = kr>D. The value of ku used by tempo2 (and 
tempo) is &d = 2.410 x 10~ 16 cm -3 pc; see Paper I for further discussion. 

Many pulsars exhibit temporal variation in the interstellar dispersion measure. tempo2 follows other software packages 
in modelling these variations by means of a Taylor series: 

DM = ^5M^(tf B -t D )", (39 ) 

n>0 

where DM'"' is the (fittable) n-th derivative of the dispersion measure and to is a user-specified epoch. 

This approach is limited to modelling variations of low to moderate complexity. At the level of accuracy demanded by 
current high precision timing campaigns, significant unmodelled variations will appear on short timescales, due to ionospheric 



5 iHellingj (1986.) quotes a figure of 36 ns, referring to the calculation of Richter & Matzncr 11983). However, that calculation incurred 
an extra pair of factors of two, one due to the fact that the calculation referred to a round trip, and one due the fact that the reflector 
was located 1 AU past the sun, rather than being a very distant source. 
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f Section I2.5.2H . interplanetary (Section 12.5,41 and interstellar electrons Foster fc Cordeslll990l) . An alternative and recom- 
mended course of action is to use the "stride fit" feature of tempo2 (Paper I), to fit a different dispersion measure to every 
group of simultaneous (or contemporaneous) multi-frequency arrival time measurements. 
tempo2 also allows for the fitting of an additional frequency-dependent delay term: 



Afdd — fcf (/ SSB ^ C 



(40) 



where k{ and £ define the scale and spectral index of the term. This term may be used in conjunction with the stride fit 
feature to model, for example, delays due to refraction in th e interstellar medium ( Foster fc Cordeslll990h or deviations from 
the cold plasma dispersion law ((Phillips fc WolszczarJll992h . 



2.6.2 Vacuum propagation delay 

The vacuum propagation delay is affected by the distance to the pulsar and the variation of this with time. The full geometric 
distance was derived in Section f2,3l The first term from equation @ to be assigned to Ais, \Ro\/c, is in fact not measurable 
to even remotely sufficient precision, but is constant and so can be dropped from the timing formula. A side effect of this is 
that epochs of measured pulse frequencies, times of periastrons, glitches and so on are in fact retarded by the initial light 
propagation time, c|iJo|. The remaining parts assigned to the interstellar delay are as follows: 

sis=fc '' + W^w- (41) 

The first term on the right hand side of equation (141 \ simply represents the displacement of the system in the initially radial 
direction. This receives potentially significant contributions from the radial veloc ity of the pulsar and the radial acceleration 
in the gravitational potential of the gl obular cluster o r Galactic environment (|Damour fc Tavloi1ll99ll) . The second term 
describes the so- called Shklovskii effect dShklovskiilll970h . whereby initially transverse motion attains a radial component due 
to the change in the direction of the line of sight caused by the transverse motion of the pulsar. That is, the pulsar moves 
tangentially while a path of constant distance describes a circle centred on the SSB. The distance between the two is the 
excess propagation length. The third term in corresponds to secular change in the size of the Shklovskii effect as the radial 
motion increases or decreases the distance, and hence changes the curvature of the line of constant distance just mentioned. 
An alternative interpretation of this term is that it is the cumulative effect upon the Shklovskii term caused by the increase or 
decrease in the apparent proper motion owing to the i ncreasing transvers e component of the initially radial motion (Section 
12.5.31 . To our knowledge, this effect was first noted bv lvan Stratenl J2003h . in the context of its manifestation as an apparent 
second derivative of the spin and orbital periods. 

These effects are often neglected in timing formulae, because as simple constant or secularly increasing Doppler shifts, 
they are inseparable from the pulse and orbital periods and their derivatives. A common procedure is to correct the apparent 
period derivatives by estimating the contributions of transverse motion and gravitational acceleration and subtracting them 
to obtain the intrinsic value fe.g. lDamour fc Tavlorlll99li ICamilo et al~lll994h . Alternatively, in certain binary systems where 
the true orbital period derivative is expected to be negligible, the measured derivative can be attributed en tirely to these 
effect s to obtain, in combination with a measurement of the transverse proper motion, an estimate of \Ro\ jBell fc Bailed 
Il99( jfl. Likewise, a measured second spin frequency derivative could, in principle, yield a measurement of the radial velocity 
l|v^aiSteatenll20n^l via the third term of equation l|41|l . 

As an alternative to modelling these effects via artificial contributions to spin and binary parameters, tempo2 can 
optionally include them directly in the following form: 

A V \\ / ,BB . \ . a ll +dshk\fJ.±f ( BB , \2 1 / , |2W,BB . \3 . 

AvP = — (t a -tposJH — (t a -ipos) + — (a M - -UHShkdotlM-Ll j (*a -tposj , (42) 

where fiishkdotj a ||> ^Shk and a M = a± -^i± are fittable parameters that specify the radial velocity, radial acceleration, Shklovskii 
distance, and component of acceleration in the direction of the transverse proper motion. The terms involving acceleration 
arise due to the expansion of k in terms of proper motion, distance and acceleration (equation|SJ. The circular dependence of 
Ayp and t^ B is resolved by iteration as described in Section f2,5.3l 

The parameters a« and dshk are inseparable, and highly covariant with a simultaneous change to the spin frequency 
derivative and orbital period derivative. It is therefore imperative that at least one of the first two parameters plus any one 
of the remaining three parameters are held fixed at some independently determined value when fitting. These parameters are 
included separately because each of them is potentially subject to external constraints. The transverse proper motion, fj,±, 
is typically immune from this covariance due to its appearance in the annual proper motion term (Section 12,5.3^ . Likewise, 
^UShkdot and are inseparable, and highly covariant with the spin frequency second derivative, so only one of these parameters 
should be allowed to vary in a fit. They are included separately because Vu may be determined from fxu (Section 12. 5. 3|l and 
one or more of the distance parameters, while a may be estimated from models of the Galactic rotation and gravitational 
potential. The first-order radial velocity term is discussed in the following section. 
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2.6.3 Transformation from SSB to BB coordinate time 

The final term contributing to Aig accounts for the special relativistic time dilation owing to the relative velocities of the 
Solar system and binary barycentres: 
,ssb _ , 

t BB -t pos = ^= tpos , (43) 

y/i- \v\y c z 

where v is the relative velocity of the frames. This transformation is applied to the arrival time of the pulse at the BB, so 
that 

tT = ta SB -Ars D -A FDD -Avp-A ES , where (44) 

2 

A ES = ^ (if SB - t pos - Aisd - Afdd - A V p) , (45) 

neglecting terms of order v 4 /c 4 and higher. Here the notation Aes refers to the fact that the term is analogous to the Solar 
system and binary Einstein delays ('Sections 12. 5. II and 12. 7. 4|l and relates to the relative secular motion of the SSB and BB. 

This transformation combines with the much larger Vu/c term of equation H'll to effect a Doppler shift between the two 
frames: 

r = ^(^ SB -W) (46) 

(47) 



y/1 -V 2 /C 2 

An unknown Doppler shift of this kind cannot be distinguished from a re-scaling of various parameters intended to refer to the 
BB frame, meaning that neither vn nor v 2 can be measured via their contribution to V. In fact, because the numerical values 
used for c and G are the same in all frames, the effect of neglecting Aes is the same as a change of unit system involving 
re-scaling of the length ma ss and time units, and ph ysical tests involving only "measured" (wrong) parameters and these 
constants will remain valid jDamour fc Deruellel [r986 ) . 

In the timing model presented here (and in contrast to that of lDamour fc DeruellellT986t) . v can in principle be deduced 
from parameters measured via other terms of equation JKJ. For this reason it may be included by specifying non-fittable 
values for the initial radial velocity vv, and transverse speed v± in the parameter file, which correct for the Doppler effect via 
equation (1 1 Lift and the following formulation of the relativistic second-order Doppler shift: 

v 2 + v 2 

A ES = - 2 J ± (if B - t poB - A ISD - A FDD - A V p) (48) 

However, it must be remembered that in practice, the uncertainty of many parameters measured in the frame of the binary 
barycentre or pulsar will be dominated by the uncertainties in vu and v±. 



2.7 The effects of a binary companion 

The effects of a binary companion are separated into several components: the geometric Roemer delay, which accounts for 
excess the vacuum light travel time due to the Euclidean displacement of the pulsar, a pseudo-delay that accounts for the 
aberration of the radio beam by binary motion, the Einstein delay (combined gravitational redshift and special relativistic 
time dilation in the pulsar frame) and Shapiro delay (gravitational time dilation in the vicinity of the companion): 

A B = A RB + A AB + A EB + A SB . (49) 

The binary model used by tempo2 i s based upon those of iBlandford fc Teukolskvl l)l976h . Da mour fc Deruelle ( 1986; 
hereafter D P) . iTavlor fc Weisber d l)l989h . Wex (1998 unpublished contribution to TEMPO; see also lLange et alJ^OOlF) and 
IWexl Jl998h . with extra terms as described bv lKopeikir] Jl995lll99rJ) . 



2.7.1 Roemer delay and Kopeikin terms 

The variation in the distance of the pulsar due to binary motion contributes four significant terms to the full geometric path 
length of equation JSJ: 

SB = h + W\ { k± b ^- r ^ b ± + ■ ( 5 °) 

The first term of equation l|5fl|l is the well-known first-order geometric delay due to the initially radial component of the 
displacement due to binary motion. The remaining terms are much smaller in magnitude and have only proven measurable 
in a few pulsars to date. Hereafter we shall refer to these collectively as the "Kopeikin terms". 
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The second term of equat ion (1501 describes change s to the apparen t viewing geometry of the orbital motion, due to the 
proper motion of the system jArzoumanian et alJll99a : lKopeikinl fl996). An approach taken in the past to account for this 
effect is to allow it to be absorbed in derivatives of the projected semi-major axis, x, and longitude of periastron, w, with 
the measured values being converted into constraints on the position angle of ascending node, il, and the inclination, i (e.g. 
ISandhu et alJll997l : iNice et al.ll200lh . Alternatively, the binary model can be mod ified to explicitly include linear changes to 
x and to as parameterised by fittable parameters Q, and i (e.s. Ivan StraterJboO^ . Here we take a third approach which we 
believe is cleaner: rather than defining the orbital parameters in a rotating reference frame, we state explicitly that the orbital 
parameters refer to the reference frame defined by the "initial" geometry at time tBB, and hence account for the effect directly 
via the second term of equation 15011 . 

The third term of equation 15011 describes the so-called annual-orbital parallax iKopeikinlll995t) . This may be interpreted 
either as the modulation of the Solar system Roemer delay (Section l2.5.3H due to the transverse orbital motion of the pulsar, 
or the modulation of the binary Roemer delay (first term of equation I5UH due to the transverse orbital motion of the Earth. 
As with the secular changes to the apparent viewing geometry, alth ough the annual-orbital parallax has been modelled in 
the past by perturbing the orbital parameters with terms involving r (Ivan Straten et alJl200ll : ISpraver et alj l2005). we choose 
maintain a consistent reference frame for the orbital parameters and instead model the effect directly via the third term of 
equation 1501 1. 

The fourth term of equation 150H is the orbital parallax iKooeikinl Il995h . which accounts for the fact that transverse 
binary motion alters the line of sight and thereby acquires an apparently radial component: this is the orbital equivalent of the 
Shklovskii effect (Section l2~6~2l . To our knowledge this is yet to be measured in any pulsar, and naturally cannot be absorbed 
in time- varying orbital parameters since the changes themselves occur on the orbital timescale. 

The contributions of the four terms of equation 1501 are included in the timing model as follows: 

Are = A RB || + A K b, where (51) 

Arbii = — and (52) 
c 

* 1 / <bb . \ , r±-b± \b±\ 2 , , 

A K B = - U a - tpos ) M-LB ■ 0± ; V -—. • (53) 

c v ' cd A op 2cd p 

where c?aop and dop are estimates of the initial distance which can be either held fixed at some independently determined 
value (or omitted from the parameter file, to neglect the effect), tied to other fittable distance parameters (e.g. d p and/or 
rfshk), or, given sufficient constraints on the relevant orbital parameters, fitted to obtain independent distance estimates. 
Likewise, /ixb can be either tied to n±, or omitted from the parameter file to neglect the effect. Note that these parameters 
in fact refer to the quantities as observed in a frame that is comoving with the BB, in contrast to corresponding distance 
and proper motion parameters in preceding sections. However, the transformation can be neglected for a fractional error of 
0(v 2 j<?) <~ 1CP 5 in those parameters, corresponding to changes in the time of emission calculation well below our 1-ns goal. 



2.7.2 Post-Newtonian orbital kinematics 



Calculation of the various orbital delays depends upon a knowledge of the displacement b of the pulsar from the binary 
system barycentre, or at least its radial component 611, depending on the demanded accuracy. Tempo2 follows the generalised 
post-Newtonian treatment of DD for the calcu lation of b, with the add ition of secular derivatives for the orbital period, 
eccentricity and projected semi-major axis after iTavlor fc WeisbereJ I 



|6| 
■To) 



A e (u) = 2 tan 



ee 



{ ei e 2 

a(l — e r cos u 
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Specifically, from equations (16-17) of DD 



where 




(54) 

(55) 
(56) 
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(59) 
(60) 

(61) 

(62) 
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k = -, (63) 
n 

e = e + e(tr-T ), (64) 

where P^o and Pi, are the initial value of the orbital period and its derivative, u>o and ui are initial longitude of periastron 
and the mean of its derivative, To is the proper time of periastron, eo and e are the initial "proper time eccentricity" and 
its derivative, and S r and 5$ parameterise relativistic deformations of the orbit. The matrices of equation 1541 in right-to-left 
order account for the inclination of the orbit to the line of sight (i), rotation about the line of sight (fi), and projection from 
a radial-transverse basis to a frame rotationally aligned to the ICRS (Section 12, 2^ , The matrices are chosen such that for 
i = 0, the angular momentum vector of the orbit is antiparallel to Ro, and fi measures the position angle of the ascending 
node with respect to e-j, in the sense of rotation into ei. The definitions of i and Q therefore correspond to astronomical 
convention if ei and ea are east and north vectors. This is the case by default in tempo2: e*i = a, = 8. Alternatively, if 
ecliptic coordinates are in use (Section I2.5.3L ei = A and = (3, so that position angles are measured counter-clockwise 
from an ecliptic meridian. 

Note that the solution to Kepler's equation 15611 involves the proper time of emission, t^ ST . This depends upon the 
equivalent coordinate time of arrival at the BB, which in turn is related to the proper time of reception at the observatory 
by all of the terms in the timing model, including those involving t^ ST . This circular dependence is resolved by starting with 
tp SI = t^ B , computing the orbital delays, making an updated estimate of tp Sr , and iterating this process until the change in 
the orbital delay is less than 100 ps. 

Using the above expression for b and writing it in terms of its transverse and radial components, the geometric delay of 
equation 15 1^ can be found explicitly. Beginning by expanding S = \b\/asin6 and C = \b\/acosd, 

S — s'mu (cosu — e r ) + cosw (l — eg) 1//2 sinu, and (65) 

C — cosoj (cosit — e r ) — sinoj (l — ea) 1//2 sin«, (66) 

where we have followed the Appendix of DD in making a minor modification to the definition of e to simplify the expressions. 

Writing equation 1541 in terms of the projections of b upon ei, and Ro as a function of S and C, 

bi = a(C sin Q + S cos f2 cos i) , (67) 

62 = a (C cos Q — S sin f2 cos i) , and (68) 

by = aSsini. (69) 
The basic Roemer delay may then be written: 

Arb|| = xS, where (70) 

x = x + x (C - To) , (71) 

and x = a/csmi is t he projected semi-major a xis as a vacuum light travel time, defined in terms of its initial value, xq (DD) 
and its derivative, x llTavlor fc Weisberslll989l) . For the Kopeikin terms, following equation 1531 

A KB = A S r + Aaop + A p, where (72) 

Asa = x (t a B — tpos) [(mi sm ^ + cos f2)C esc i + (^11 cos SI — ^12 sin D.)S cot i] , (73) 
x 

Aaop — ~ -3 [(r • ei sin O + r ■ £2 cos Q)C esc i + (r ■ e± cos fl — r ■ £"2 sin Q)S cot i] , and (74) 

oaop 
2 

Aop = ^— { C 2 esc 2 i + S 2 cot 2 i) . (75) 
2o!op 

The proper motion components are defined in terms of fittable parameters that can optionally float independently to the 
annual proper motion parameters (Section 12. 5. 311 . in order to allow the latter to be separated from the secular change of 
binary viewing geometry. For equatorial coordinates, hi — /i a B and 112 = /Mb, while for ecliptic coordinates /ii = /j,\b and 
M2 = M/3B- In agreement with astronomical convention, the parameter Q corresponds to the position angle of periastron in 
terms of eastward rotation from north, while the parameter i is defin ed such that i = cor responds to a binary orbital angular 
momentum vector dir ected towards the observer. As pointed out bv lSolaver et alJ l)2005t) . these differ from the definitions of 
iKopeikinl Jl995t [19961^ . We note that some caution is required in this area owi ng to variations in d efinit io ns in the literature: 
for ot h er examples of use of the unconventional definition of i see e.g. DD; iDamour fc Tavlorl Jl992l) : IWeisberg fc Tavloil 
J2002h : lstairs et al l J2004h . 



2.7.3 Aberration 



Owing to the relative transverse velocity of the pulsar and the Earth, the direction of the observer as seen from the pulsar 
differs from that which would be measured in a frame co-moving with the observer. The direction is aberrated according 
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to the Lorentz transformation that relates the two frames. Under the assumption that the source of pulses is the periodic 
rotation of an emission beam, any such change in the direction of the observer alters the rotational phase corresponding to 
radiation received at a given time. Following DD (equation 27), in TEMP02 this is converted to an equivalent time delay: 

Aab = A {sm[uj + A c (u)] + e sin a;} + B {cos[ui + A c (u)] + e cos lu} , (76) 

where A and B axe parameters related to the orientation of the pulsar spin axis and the size of the orbit (see DD equations 
38-39 and DD section 3.2). As discussed by DD, these parameters are highly covariant with the other binary parameters, and 
may only be separated on the time scale of geodetic precession. Note, here the relative motion of the observatory and the BB 
is neglected. 



2.7.4 Post-Newtonian delays 

In addition to the modified Roemer delay of equation 151H and the aberration pseudo-delay, under the post-Keplerian formalism 
of DD there are two additional delay terms: the Einstein delay and the Shapiro delay. 

The Einstein delay is the difference between the proper time of emission and the coordinate time in the quasi-inertial 
frame of the binary barycentre: 

t BB =t pBI + A EB - (77) 

This is due to the combined effect of gravitational redshift and time dilation. Following DD and in contrast to the tempo2 
treatment of the Solar system Einstein delay, the binary Einstein delay specifically excludes a linear term by scaling t pBT in 
such a way that the orbital period is numerically the same in either time scale. It is expressed in a theory-independent manner 
in terms of a timing model parameter, 7; from equation (19) of DD: 

Aeb = 7 sin it. (78) 

The Shapiro delay is caused by the curvature of space-time in the gravitational field of the binary companion. After 
equation (26) of DD: 

.1/2 



A SB = -2rlog \ 1 



sin a; (cosu — e) + (l — e 2 ) 



cos ui sin u 



(79) 



Here the fittable parameters are the theory-independent "range", r, and "shape", s. Under general relativity, s — sini and 
r = Giri2/c 3 , where 7712 is the mass of the binary companion. In tempo2, r is expressed in units of Tq — GMq/c 3 , so that its 
value is numerically equal to the GR companion mass in units solar masses. (Tq is half the light travel time across the Solar 
Schwarzschild radius.) An alternative parameterisation is available, where the parameter s is replaced by the "shapmax" 
parameter, z s = — ln(l — s), This derives fro m a modification to tempo designed to avoid difficulties in fitting highly inclined 
orbits near | sin i \ — 1 l|Kramer et alJl2006ah . In this case, the above relation applies after substituting s — 1 — exp( — z s ). 



2.7.5 Assuming general relativity 

The post-Keplerian parameters of Sections 12.7, 2l and l2, 7.41 f A. Co, 5 r , 5g, 7, s, 7712) can either be specified directly, or they can 
be derived from a reduced set of physical parameters based on the assumption that the predictions made by general relativity 
are correct. These parameters are and M = mi + m,2, the total system mass (in addition to the Keplerian parameters, 
Pbo, <^o, e, To, and xq). After equations 8.48-8.55 of Lorimer & Kramer (2005), 



• < ;<; 0^2/3 5/3 M 2 ^ 3 

u = 3Tq n — , (80) 



7™ - r^n-^e — — — , (81) 



.gr _ T 2/3^-i/3„ra2 (mi +2m 2 ) 
r GR = T Q m 2 , (82) 

GR • • TH-1/3 2/3 M 2 / 3 , 

s — smt — T a ' n ' x , (83) 

w m2 

cGR _ ^2/3 2/3 3m 2 + 6mim 2 + 2m^ 

°r - 1 Q n M 4/3 ' 

cGR ^2/3 2/3 l m i + 6mim 2 + 2ml , , a _. 
Se = T e ' n ^ , and (85) 

• GR 192tt 5/3 5/3 m im2 1 + 73e 2 /24 + 37e 4 /96 

P » - ^ T n MTTF (1 _ e 2)7/2 ■ ( 86 ) 

Although these equations account fully for the general relativistic rates of change of orbital period and longitude of 
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periastron, extra variation in these parameters is allowed for via the parameters P bx and lu x 



n GR = ^- + ^^ p 7 - " and (87) 

no p b0 

cu GK = u; +(k Gn +^A e (u). (88) 

Note that ui x is not intended to model the secular change in orbital viewing geometry due to proper motion: this is modeled 
by a separate term in the timing model, causing u) to be defined in a consistent, non-rotating reference frame (Section 12.7, II . 

Under the condensed parameterisation of this section, the orientation of the spin axis may be separated from the ori- 
entation of the orbit normal vector to yield a physical parameterisation of the aberration delay (DD) . Following equations 
8.56-8.57 of Lorimer & Kramer (2005): 



2n n(Pt K + P bx )(tr -T ) 



A ° R = TTvm TR — sin X esc ^, (89) 

(2tt) 2 / 3 pPl'\l - e 2)l/2 M 2/3 



-1/3 



T 

B GK = - _ ., ,, cos i cos x esc £, (90) 

(27r) 2 / 3 vPl /3 {l - e 2 )i/2M2/3 

where v is the pulsar spin frequency (Section 12.81 . £ is the inclination of the pulsar spin axis, and % i s its position angle 6 , 
relative to the position angle of ascending node, Q. Here \ and csc£ are the fi ttable parameters. Note that because cosi is 
only ambigiously constrained by the Shapiro delay term (cosi = ±yl — sin 2 i), \ too is subject to ambiguity. To evaluate 
equation 1901 . tempo2 assumes that cosi > 0. Values of \ determined on this basis are therefore indistinguishable under 
the transformation \ ~~ * n — X- Only if i is measurable via the Kopeikin terms can this degeneracy be broken. A second, 
unbreakable degeneracy exists due to the fact that aberration is only sensitive to the projection of the transverse orbital 
velocity upon the spin axis: £ — » n — f . 



-1/3 



2.7.6 Nearly circular orbits 

For low-eccentricity orbits, uo and To are highly covariant, which complicates the analysis of the uncertainty in these parameters. 
This problem can be avoided by parameterising the motion in terms of the Laplace-Lagrange parameters, 

r/ = esincj and (91) 
K = ecosoj, (92) 

and the time of ascending node, 

T asc ee To - u/n (93) 

(Wex 1998 unpublished contribution to tempo as "ELL1" binary model). Following lLange et, alJ l|200lfl . the following ap- 
proximation to equation l|7()|l applies, to order e: 



LL nLL 



i fiS|| 



X 



S , where (94) 

5* LL = sin $ + | sin 2$ - ~ cos 2$, (95) 

$ = n(tr- Tasc), (96) 

X LL = X0 + X (if" - ^asc) , (97) 

V = T] + i) {C 1 - T asc ) , and (98) 

K = Ko + « (C - Tasc) ■ (99) 

Here 770, ??, Ko, k and T asc replace eo, e, u>o, a; and Tp as fi t table parameters, while 5e and S r are neglected. 

The original ELL1 model as presented bv lLange et alJ i200ll) did not include the Kopeikin terms (Section l2.7.1t . In order 
to include them, an expression for C (equation 1661 is required. To the same accuracy as equation 1951 . we find 

C LL = cos$ + ^cos2$+ I sin 2$, (100) 

neglecting a constant offset of 3aAc/2 (analogous to an offset of 3xr//2 dropped by lLange et allbOQlh . Substitution into 
equations i|73l , 174^ and l|75|l yields the Kopeikin terms under this parameterisation. 

Under the Laplace-Lagrange parameterisation, the Shapiro delay is computed as follows l)Lanee et alJl20o"l|) : 



6 These are related to the rj and A of lDamour fc Tavloil ll992t) and Lorimer & Kramer (2005) via \ = ~V> ? = 7T — A. The difference is 
to ensure the use of the conventional definition of position angle and inclination, see Section I2T31 
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Asb = -^^ln(l-sinisin$), (101) 

c 3 

where either of the parameterisations (s, z s ) of sin i of Section \2. 7, 41 may be used. 
The aberration delay becomes (Wex 1998 unpublished): 

Aab = .4 sin* + Bcos$. (102) 

The Einstein delay is neglected if this parameterisation is chosen. 

2. 7. 7 Main sequence companions 

The orbit s of pulsars with main sequence binary companions show deviations from a Keplerian orbit attributed to spin-orbit 
coupling jLai et aJll995h . tempo includes two binary models for taking these into account: the BTJ model, which induces 
steps ("jumps") in several parameters at specified epoc hs (typicall y, at each periastron), and the MSS model, which includes 
orbital and secular variations to the inclination angle (IWexl 1998). The insertion of jumps is straightforward and need not 
be discussed further. The MSS model is included in the tempo2 model by changing the manner in which x alters x, and 
including second period derivatives in x and ui: 

«"» = , + ^ + ^ PSr - r °> 2 ,and (103) 
n 2 

v MSS = u + kA e (u) + *( tPST - To) \ (104) 

2.7.8 Emulating other models 

As noted in Section 12.7.21 the orbital kinematics represent an inverse problem, because the motion is parameterised by the 
pulsar proper time, which itself depends on the orbital propagation delay. The standard approach in TEMPQ2 is to itera te until 
convergence. However, in order to reproduce results from the implementations of the IfBlandford fc Teukolskvl <ll976l) or DD 
models, a facility is provided to alter this behaviour. Both of these models begin with the zeroth-order approximation that 
tg Sr = t BB to compute the various delay terms, then add correction terms to the zeroth-order approximation to the combined 
Roemer and Einstein delay, A RE , which is expressed (after DD equations 46-48) as follows: 

A RE = ob(cos U — e r ) + (f3 + 7) sin U, where (105) 

qb = sinoj, (106) 

c 

a sin i / 2 \i/2 , , 1f w. 

f) — (1 — eg) cos to , and (107) 

c v ' 

U-e T sinU = n(t BB -T ). (108) 

Note that the Kopeikin terms (equation I53H are neglected. 
The expansion provided by DD is 

2 . 1 „ 2 . . n'" 1 e sin U 



a dd A o 

fl RE — ^RE 



nA RE + n A RE + -n A RE A RE - — n A RE A RE 

V / 2 2 1 — e cos U 



1 — e cos U ' 



, where (109) 

(110) 



Are = -Qb sin U + (0 + 7) cos U, and (111) 

A RE = -qb cos U - (p + 7) sin U. (112) 

The approximation provided bv lBlandford fc Teukolskvl lll976l) is equivalent to: 

All = AR E (l-nAR E ). (113) 

In addition. iBlandford fc Teukolskvl (1 19 7&i neglected Abs, Aba, 5 r and 8$, and the orbital modulation of u: that is, in contrast 
to equation 1621 . 

Lu +u(t psl -To). (114) 



Keeping the DD formulation of ui but the BT formulation of the other terms allows for the emulation of the BT+ model of 
tempo foamour fc Tavlorlll992h . 

Additionally, it is possible to emulate the tempo implementation of the Laplace-Lagrange parameterisation (Section 
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I2.7.6L which made an approximation to the Roemer delay analogous to the DD treatment: 



A 



i aBLLI-O' . ( aELLI-O'A" 3 . 1 2 a ELL1-0 a ELL1-0" 

1 - nA RB + (^nA RB j + -n A RB A 



where 



sin $ + — sin 2$' 



o ?? 



cos 2$ 



a ELL1 — 0' A 

A RB = a;cos$, 



ELL1-0" 
RB 



= — csin$, and 
$ = n(^ B -T asc ). 



(115) 

(116) 

(117) 
(118) 
(119) 



2.8 Pulse phase model 

The final component of the timing model is the evolution of the phase of the pulse sequence, relative to the proper time i psr 
of the pulsar centre of mass (equation |7J. In most cases, a simple Taylor series expansion is used: 

L — ' nl 

n>l 

The frequency derivative terms (u^ n ') are the fittable parameters, while the epoch, tp, at which <j> = v is set by the user. 
Absolute phase alignment is achieved by means of 4>o, which is not specified directly. Instead, in order to protect against 
changes in tp or the choice of barycentric time scale, <j)o is defined in terms of a reference time of arrival for a specified 
observing site and observing frequency. This time of arrival is chosen such that it is close to the centre of the time span of 
the fitted data, and has a predicted pulse phase that is exactly integer. 

Note that under the assumption that the basic pulse train is a result of a rotating emission beam, the time of zero phase is 
in fact dependent on the position of the observer. Of all the sources of relative motion (Earth orbital, binary orbital and secular 
motion and acceleration), only the secular motion of the pulsar (or binary barycentre) is significant given our accuracy goal 
of 1 ns. However, since the transverse speed is nearly constant, so is the rate of change of apparent pulse phase, meaning that 
the effect cannot be separated from the basic pulse frequency, v. For this reason, the effect is neglected in tempo2. Likewise, 
the emitted beam is subject to aberration due to the relative velocity of the pulsar, however the constant component of this 
velocity contributes only a constant phase shift that is absorbed in <f>o- On the other hand, the binary motion induces a 
variable aberration, which for simplicity is included in As as a pseudo-propagation delay Aa (Section l2.7.1H . 

Though the pulse frequency of most pulsars shows a continual gradual slow-down, attributed to magnetic braking of a 
spinning neutron star, certain pulsars exhibit anomalies known as "glitches". These events are generally well- modelled by 
permanent increments to the phase, frequency and frequency first derivative, in addition to a frequency increment that later 
decays exponentially to zero: 



A0 + Au(t^ si - y + |A;>(tp sr - t g f + fi - e -(*r-tg)/^ Av t {tr r - t e ) i<r > t s 
o tr r < u. 



(121) 



An arbitrary number of glitches may be modelled in this manner, the corresponding values of <j)g being added to the basic 
Taylor series of equation l|12()|l to give the predicted evolution of pulse phase with time. 



3 ACCURACY ESTIMATES 
3.1 Geometric propagation delay 

Errors are introduced to the geometric part of the propagation delay in two ways: through the neglect of various terms of 
equation @, and through approximations made in the computation of other terms. Cases of the latter are discussed below. 
The neglected terms are the O(|i?o| -3 ) terms, and all those involving the expansion of the third, second, and fourth terms 
in the second set of parentheses, except for the secular and annual changes to the Shklovskii effect and the secular change to 
the annual proper motion. Based upon measured parameters and a timing campaign spanning 20 years, we have estimated 
the likely magnitude of the neglected terms for all pulsars in the current ATNF pulsar catalogue having spin periods shorter 
than 30 ms (slower pulsars being unlikely to yield high precision timing measurements). The largest of the neglected terms 
are the change in annual proper motion due to transverse acceleration (~ 1 ns for PSR B1620— 26), second-order secular 
changes to the apparent binary orbital viewing geometry (~ 600 ps for PSR J0437-4715), the Od-RoT 3 ) terms involving 
secular motion (~ 600 ps for PSR B1257+12), the binary modulation of the Shklovskii effect (~ 300 ps for PSR J0437-4715) 
and the second-order annual proper motion (~ 100 ps for PSR J0437— 4715). 

It is possible that pulsars yet to be discovered may yield larger errors than the examples listed here, if they are in wider 
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Figure 1. Round-trip time transfer error involving independent clock monitoring systems at Parkes observatory. The terms in each 
correction are UTC(GPS)-UTC(PKS) ("Totally Accurate Clock" system), UTC-UTC(GPS) (Circular T), UTC-UTC(AUS) (Circular 
T) and UTC(AUS)-UTC(PKS) (GPS common view system). 

binary orbits, are closer to the Earth, have larger proper motions, and/or experience greater secular acceleration. The current 
sample is not likely to be biased in relation to proper motion, and any bias in relation to distance will tend to select closer 
pulsars due to the dependence of detection upon flux density. Highly accelerated pulsars are difficult to detect, however the 
level of acceleration at which this becomes problematic is many orders of magnitude larger than may be experienced a the 
Galactic or intracluster environment. Binary acceleration also hampers detection, but for the known population the lack of 
any strong dependence of companion mass upon orbital separation means that this selection effect biases the sample towards 
large orbital separations rather than small. However, it is not inconceivable that the present sample has missed an entire 
class of millisecond pulsars that experience strong binary acceleration due to very massive companions in wide orbits. The 
discovery of such pulsars may necessitate the re-inclusion of some of the dropped terms of the geometric propagation delay. 

3.2 Clock Corrections 

Clock corrections applied by tempo2 may suffer from time-transfer errors present in the clock offset files provided by the 
user, owing either to incorrect tabulated values, or to variations in the clock offset on time scales shorter than the tabulation 
interval, over which linear interpolation is performed. In most instances, the dominant source of these errors will be the 
first step of correction, involving the local observatory clock, since observatory time standards and time transfer systems are 
typically less accurate than those involved in GPS and UTC time keeping. As an example of the level at which time transfer 
errors may occur, we note that differencing of two independent time transfer systems at the Parkes observatory leads to an 
error term with a root-mean-squared (rms) deviation of 8 ns over the last three years (Fig. [I). Time transfer to UTC from 
GPS or one of the nation al time standards contributing to UTC also leads significant error, with a typical rms of ~ 5 ns 
llLewandowski et al.ll2005l) . 

In addition to clock correction errors, the "corrected" time of arrival, as measured against a specified realisation of TT will 
differ from that measured against an ideal TT, owing to instability in the contributing atomic clock systems. An indication of 
the level at which instability contributes can be obtained by means of the difference between the latest revision of the BIPM 
realisation of TT, TT(BIPM2005), with the 2001 version of the same timescale. As shown in Fig. [5] differences exist at the 
~ 10 ns level, even after allowing for the absorption of a fraction of the variation by other terms in the timing model. 

3.3 Solar system Einstein delay 

The time dilation integral of llrwin fc Fukushimal \ 19991) (equation includes only terms of order 1/c 2 , plus a mean rate 
correction for higher-order terms. Inaccuracies in the Solar system ephemeris lead to an estimated uncertainty of 100 ps in the 
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Figure 2. Difference between two realisations of TT, TT(BIPM2005) and TT(BIPM2001) before (solid line) and after (dashed line) 
removing a fitted function including a quadratic and an annual sinusoid, to account for the presence of fitted functions of this form in 
the timing model. The plot spans the final five years covered by TT(BIPM2001). 



time integral, and 10 -16 in its derivative. The remaining higher-order terms have a maximum amplitude of 33 ps iFukushimal 
Il995l) . The effect upon the spatial coordinates of the 4-dimensional transformation relating the GCRS and BCRS is neglected. 
The neglected length contraction is of the order 10 -8 s/c, corresponding to up to 200 ps of excess propagation time. 



3.4 Atmospheric delays 

Equation 1111 predicts the hydrostatic zenith delay to an accuracy of 0.5 mm under conditio ns of hydrostatic equilibrium, 
with variations of up to 20 mm possible under extreme weather conditions foavis et alJll985l) . These correspo nd to timin g 
errors of 1.7 - 67 ps • esc 0. In addition, the accuracy Niell mapping function degrades at small elevation ang les jNiel]|ll996h . 
leading to ~ 30 ps of error at = 5°. Should the user not provide tabulated surface atmospheric pressure, the use of a 
constant value leads to errors of the order of 1 ns ■ esc 0. Likewise, lack of tabulated values for the zenith wet delay results in 
error at the level of ~ 1 ns- cscO. 

As noted in Section f2.5.2l the dispersive component of the ionosphere causes propagation delays of the order of 
7—130 ns -(//l GHz) -2 , where / is the observing frequency. Such effects must be removed by fitting for a time- variable 
dispersion measure with the aid of multi-frequency observations. 



3.5 Solar system Roemer delay 

The Roemer delay depends on the position of the observatory in the terrestrial frame, provided as input to tempo2. An 
accuracy of 30 cm is needed to meet the 1 ns criterion. In many cases the position is known to an accuracy of a few 
centimetres, by means of VLBI or other geodesy techniques. In addition, most sites are subject to continental drift at a rate 
of a few centimetres per year, which is neglected by tempo2. This may be remedied in future when the sensitivity and time 
span of observational data dictates it. 

The transformation of the observatory position from the terrestrial to celestial reference frame depends on a knowledge of 
the Earth orientation parameters. Our requirement of 1 ns precision corresponds to maximum 10 mas of angular displacement 
on the surface of the Earth, or 600 /is in UT1, the time scale tied to Earth rotation 7 . tempo2 uses the "C04" Earth orientation 
parameter series of the IERS, which exceeds the required precision for all dates after 1980. The effects of diurnal and subdiurnal 



Errors in UT1 relate in fact to a position offset, and should not be confused with errors that appear directly in the timing model (such 
as clock corrections). This is why the precision demanded of UT1 is some 600,000 times higher than the required timing accuracy. 
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tides are removed from the series, resulting in an additional error of at most 0.7 mas in the location of the pole, and _|Jj fis 
in UT1. The IAU2000B precession-nutation model used by tempo2 is accurate to 1 mas. 

As noted in Section 153)21 tempo2 applies transformations to the input and output values of the Solar system ephemeris 
that are specific to recent JPL ephemerides. Should a different ephemeris be specified, it is important that the offset of 
the input time parameter , and the scale of input time and output spatial vectors match those of DE405, as specified by 
llrwin fc Fukushimal l|l99flF) . Likewise, it is important that the ephemeris coordinate frame is rotationally aligned with the 
ICRS, for while the gross effect of a frame rotation on r® ■ Ro could be absorbed in the fitted pulsar position, the resultant 
changes in Ro will induce an error in s ■ Ro wit h a diurnal signatur e. For example, the coordinate frame of the earlier DE200 
ephemeris is offset from the ICRS by ~14 mas jFolkner et al.lll994h . resulting in up to ±1.4 ns of error. 

Ultimately the accuracy of the Roemer delay hinges upon the chosen Solar system ephemeris. Preliminary investigations 
indicate that errors in the Solar system ephemeris may be significant and indeed measurable in pulsar timing data. 



3.6 Solar system Shapiro delays 

tempo2 includes the first-order Shapiro delay due all bodies for which the excess delay exceeds 0.5 ns (see Paper I). The next 
largest effect is due to Mars, with 60 ps of delay. The second order term is included only for the Sun. The next largest effect is 
due t o Jupiter, with a maximum amplitude of 5 ps. Further terms appear for the Sun at the 500 ps level jRichter fc Matznerl 




3.7 Interstellar propagation delays 

The treatment of tempo2 is precise in the limit of cold plasma dispersion with negligible refraction. In practise, turbulence 
in the interstellar medium leads to stochastic refraction and scattering of incident radio waves, giving rise not only to familiar 
scintillation effects, but also to frequency-dependent variable propagation delays llFoster fc Cordeslll990l) . Such delays arise 
both due to differences in the propagation path length, and also due to induced error s in the computed Roem er delay, owing 
to the variations in apparent source direction. According to simulations conducted bv lFoster fc Cordesl l)l990|l . the amplitude 
of these terms in the selected pulsars is of the order of hundreds of nanoseconds, for observations conducted at A ~ 20 cm. 
For future, more sensitive telescopes, these errors will almost certainly dominate over radiometer noise unless steps are taken 
to correct for the effect or avoid it by observing at higher frequencies. 



3.8 Transformation from SSB to BB frame 

The transformation from the relativistic coordinate frames of the SSB and BB is modelled as a rotationless Lorentz trans- 
formation. This is a pure special relativistic transformation that assumes that the space-time metric tensors are completely 
determined by bodies in the Solar system and binary system respectively. Under general relativity this assumption will be 
violated by the presence of gravitational waves from external sources, causing unmodelled effects in the pulse times of arrival. 
Indeed, the detect ion of the gravitational wave background is one of the main aims of contemporary pulsar timing campaigns. 
Ijenet et alJ J2005h conducted simulations of a number of different hypothetical datasets to test the detectability of a gravita- 
tional wave background at the level predicted by theory. A significant detection was predicted for several simulated datasets 
of a standard optimistically within reach of current observing programs, for which the minimum detectable systematic effect 
(as indicated by the rms measurement error divided by the square root of the number of independent observations) is one to 
several nanoseconds. 



3.9 Binary orbital effects 

The principal source of error in the formulation of the timing delay due to the presence of a b inary companion is the use of 
equations of motion derived at only the first post-Newtonian order l|Damour fc DeruelleHlflSfih . To date no complete timing 
formula ac curate at th e second post-Newtonian (2PN) order has been constructed, although the orbital motion has been 
derived bv IWexl Jl995h . including a proper-time representation that thereby provides the 2PN Einstein delay. The largest 
discrepancies are due to the 2PN contribution to the rate of periastron advance, and the effect of spin-orbit interaction on 
the periastron advance and orbital orientation JPamour fc Tavlorlll992l) . though these effects can be absorbed i n changes t o 
u) and x. The 2PN corrections to the quasi-periodic Einstein delay are of the order 80 ns for PSR B1913+16 l|Wexlll995tl . 
Additional unmodelled effects include the modification to the Shapiro delay due to the gravitomagnetic field of a spinning 
companion JLaeuna fc Wolszczanlll997t). and the effect of gr avitational light bending on the apparen t rotational phase of th e 
pulsar beam iSchneicterl^90nDoroshenko fc Kopeikmlll995l) . which cannot be separately measured dWex fc Kopeikmlll999l) . 
For double-neutron-star binaries, all of these effects are quite large compared to our 1 ns accuracy goal, especially for nearly 
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Figure 3. Difference between timing models of TEMP02 and TEMPO owing to the use of different numerical approximations to the Solar 
system time dilation integral. 

edge-on orbits, however it so happens that the binary evolutionary history of such systems rarely results in a pulsar that is 
spinning rapidly enough to provide timing measurements of high accuracy. 

Only since th e discovery of PSR J0737— 3039A, a 22-ms pulsar in a highly inclined compact orbit with a second pulsar 
jLvne et, al.ll20oJl . has the measurability of these effects been seriously considered. It is expected that the 2PN contribution 
to Co will be measurable in coming years JrCramer et alJ l2006b1 . Although this effect is easily accommodated in the existing 
binary model, we understand that efforts are underway to develop a fully consistent 2PN timing formula. Incorporation of this 
formula in tempo2 would be highly desirable for assurance that the model will continue to be accurate at the level demanded 
by ever-decreasing measurement errors and potentially more extreme binaries discovered in coming decades. 

It should be noted that owing to approximations made in the orbital kinematics, the La place-Lagrange parameterisation 
is only valid for orbits of low eccentricity. The error term is of order xe llLange et alJbOQll) . 



4 DIFFERENCES FROM TEMPO 
4.1 Tropospheric delays 

tempo makes no attempt to model non-dispersive atmospheric delays. These are dominated by the hydrostatic delay (Section 
12,5.20 . which is a mainly diurnal term of up to tens of nanoseconds in amplitude, depending on the elevation angle distribution. 



4.2 Einstein delay 

To com pute a version of the time dilation integr al of equation tempo uses the series appr oximation of lFairhead fc Bretagnonl 
Jl99(J) . since shown to be significantly in error jFukushimalll995l : llrwin fc Fukushimall 999 ) . An example of the discrepancy is 



shown in Fig. ED which plots the difference in barycentric time as computed by tempo2 and tempo. The comparison is made 
in terms of the now-obselete Barycentric Dynamical Time (TDB, see below), as used by tempo. The difference consists mainly 
of a constant offset of ~ 60 ns, an annual term of ~ 2 ns amplitude, and a small noise-like contribution due to truncation of the 
internal Chebyshev representation used by tempo. The constan t offset is immaterial since the zero point of TDB has never 
been well-defined. It arises for the most part due to the use by llrwin fc Fukushimal Jl999ll of an err oneous estimate of the 
offset between TDB and T cp h, with an additional ~ 6 ns due to software errors in early versions of the lFairhead fc Bretagnonl 
l|l990ll Fortran code, which only manifest on certain computer system s (P. Wallace, personal communication). The annual 
term, probably due to inaccuracies in the lFairhead fc Bretagnonl l)l990f) series, will act to corrupt the source position at the 
0.5 microarcsec level, as fitted by tempo. 

It should also be noted that the time integral actually computed bv lFairhead fc Bretagnonl Jl990l) and llrwin fc Fukushimal 
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( 1999) includes a negative term in the integrand designed to give a long-term zero mean value of the integral. Unlike tempo2, 
tempo does not add a counteracting linear term to the result, causing an overall difference in the scale of the barycentric 
coordinate timescale. In the case of tempo the resultant time refers to the T ep h timescale (Sect ion 12. 5. H . which is approxi- 
mately equivalent to the ill-defined IAU predecessor to TCB, known as TDB (see lstandishlll998h . The continued use of T ep h 
or TDB as a barycentric coordinate time for pulsar timing is deprecated, since under this condition the scale of units of fitted 
parameters involving length and/or time differ subtly from the standard SI definition and fail to comply with IAU Resolutions: 
see Paper I. 



4.3 Earth orientation 

In its calculation of the observatory positio n in the barycentric frame, tempo uses the IAU 1976 precession llLieske et alll977h 
and IAU 1980 nutation JSeidelmanr]|l982h theories, which are in error at the 50 mas level. More significantly, it also neglects 
polar motion, which amounts to up to ±300 mas, corresponding to up to ±35 ns of timing error, with a variable diurnal 
signature. The polar motion error is dominated by two periodicities, one with a period of one year (due to seasonal influences 
on oceans and the atmosphere), and the other with a period of 435 days (the Chandler wobble, due to free precession)- see 
Paper I. Depending on the site position, source direction and hour angle distribution, part of the annual polar motion may 
enter the timing model as an annual term at the level of up to ~ 10 ns, corrupting the fitted pulsar position at the level of a 
few microarcsec. 



4.4 Shapiro delays 

tempo includes only the first-order Shapiro delay term due to the Sun. Terms of up to 180 ns for the planetary Shapiro 
delays, and 9 ns for the second-order term of the Solar Shapiro delay were omitted. 



4.5 Dispersion 

tempo includes some small errors in the transformation of the observing frequency to the barycentric reference frame. Firstly, 
the computation of the site velocity (equation I27|l neglected polar motion, precession and nutation. Polar motion changes 
the direction of the velocity vector x by up to 300 mas, altering its projection on the line of sight by up to ~ 10~ 10 c and 
corrupting the dispersion term Ars of equation 13711 by a fractional amount ~ 10 -9 ' 5 . More significantly, precession moves the 
vector by ~20 arcsec yr _1 , leading to a fractional error in Ais with an amplitude increasing at a rate of up to ~ 2 x 10~ 8 yr _1 . 
This will manifest as an apparent annual variation in the dispersion measure, with an amplitude equal to ~ DM • 10 -8 t/yr, 
where t is the time elapsed since Julian year 2000.0. For single-frequency datasets, the effect would be absorbed by the annual 
proper motion term, corrupting it at the ~ 10/ias yr -1 level. 

Secondly, tempo neglected to account for the Einstein rate (Section l2.6.1|l . The use by tempo of TDB instead of TCG 
(Section 14. 21 removes a large positive mean from the Einstein rate, leaving residual rate errors at the 10 -9 ' 5 level, mostly 
modulated with an annual periodicity. This corresponds to fractional errors in Ais at the 10 -9 level. 



4.6 Secular and binary motion 

The full geometric time delay of equation © contains several terms involving the displacement of the pulsar due to secular 
motion (fc). Of these, tempo includes only those terms related both to secular and annual motion. In doing so, it neglected 
direct Doppler shift (which includes both classical and relativistic terms) due to radial velocity and its change in time due to 
radial acceleration, and the Shklovskii effect and its change over time due to the radial motion and transverse acceleration. 
None of these omissions is necessarily problematic, because all can be absorbed in a redefinition of other measured parameters. 
In tempo2 they are provided directly so that advantage may be taken of additional constraints that may be available: for 
example, if the intrinsic spin orbital period derivative is known to be small, it can be held fixed at zero and the distance 
determined directly via the Shklovskii effect. 

The public distribution of tempo also omits the Kopeikin terms, however the secular changes to the viewing geometry 
may be absorbed in other parameters, while modified versions of tempo exist that include the annual-orbital parallax. The 
orbital parallax, not modelled in any known version of tempo, contributes a timing term of up to 30 ns in amplitude. 



5 SUMMARY AND CONCLUSIONS 



We have presented a model for pulsar pulse times of arrival that exceeds all previous efforts in its accuracy. This timing 
model is the basis of tempo2, a new software package for pulsar timing. The goal of tempo2 to provide accuracy at the 
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l-nanosecond level is largely satisfied by the model presented here, the exception being pulsars that are part of extremely 
relativistic binary systems. Further theoretical work on second- and higher-order relativistic effects will be required before 
progress can be made in this area. Fortunately, the pulsars for which the most precise measurements are possible tend not to 
be part of such relativistic binaries. 

Terrestrial clock instability, Solar system ephemeris errors and the local gravitational wave background stand as the likely 
sources of significant systematic error. Pulsar timing array campaigns should eventually be able to provide measurements of 
these effects, allowing for their removal from the residuals and more importantly providing a valuable insight into the causative 
physical processes. 

With the next generation of radio telescopes, such as the Square Kilometre Array, achievable measurement error will be 
reduced in some cases to below 10 ns. For systematic effects, the timing model presented here should remain applicable in most 
cases, and will present a useful starting point for extension if necessary. However, at this level of precision previously negligible 
stochastic effects will come into play. Possibly the most problematic of these will be the effect of interstellar scattering. 
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APPENDIX A: TABLES OF PARAMETERS AND VARIABLES 
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Table Al. Parameters of the timing model 

Name Description Section 

Algebraic ASCII 



a 


R A 


Right ascension of pulsar (ICRS) 


12.5.31 


P 


ELAT 


bcliptic latitude ol pulsar 


12.5.31 


7 


GAMMA 


Amplitude ol binary tmstein delay 


12.7.41 


A(p n 


GLPH_n 


Pulse phase increment for glitch n 


12.81 


Au„ 


GLF0_n 


Pulse frequency increment for glitch n 


12.81 


Au„ 


GLFl_n 


Pulse irequency derivative increment lor glitch n 


12.81 


Al/fn 


Glf0d_7i 


Amplitude ol decaying pulse Irequency term lor glitch n 


12.81 


6 


DEC 


Declination ol pulsar (ICRS) 


12.5.31 


&e 


DTHETA 


Rclativistic longitudinal deformation ol orbit 


12.7.21 


S r 


DR 


Helativistic radial deformation of orbit 


12.7.21 


eo 


ECLOBL 


Mean obliquity of the ecliptic at J2000.0 


12.5.31 


c 


FDDI 


Exponent of frequency-dependent delay 


12.61 


no 


EPSl 


First Laplace- Lagrange parameter at i psr = To 


12.7.61 


V 


EPSlDOT 


Time derivative of first Lap lace- Lagrange parameter 


12.7.61 


Kq 


EPS2 


Second Lap lace- Lagrange parameter at £ psr = Tq 


12.7.61 


k 


EPS2DOT 


Time derivative of second Lap lace- Lagrange parameter 


12.7.61 


A 


ELONG 


Ecliptic longitude of pulsar 


12.5.31 


/-< a 


PMRA 


Proper motion in right ascension 


12.5.31 


MaB 


KPMRA 


Proper motion in right ascension for Kopeikin term 


12.7.21 


M/3 


PMELAT 


Proper motion in ecliptic latitude 


12.5.31 


M/3B 


KPMELAT 


T) J_" " l - i ■ l , ■ i l r t/ " 1 ■ j_ 

Proper motion in ecliptic latitude lor Kopcikm term 


12.7.21 


W 


PMDEC 


Proper motion in declination 


Iz.o.ol 


M5B 


KPMDEC 


Proper motion in declination for Kopeikin term 


I2T21 


MA 


PMELONG 


Proper motion in ecliptic longitude 


123^1 




KPMELONG 


Proper motion in ecliptic longitude for Kopeikin term 


I2T21 


Mil 


PMRV 


Radial proper motion 


123^1 


„(n) 


F/l 


n-th time derivative of pulse frequency at t psr = tf req 


EH 


CSC £ 


SPINCSCI 


Cosecant of inclination angle of spin axis 


|2~731 


n 


PX 


Annual parallax 


l2~Ql 




GLTD.n 


Time constant for glitch n 


nu 


00 




Pulse phase at t pBT = if roq (implicit floating parameter) 


l2~8l 


X 


SPINPA 


Position angle of spin axis w.r.t Q 


|2~731 


n 


KOM 


Position angle of ascending node 


12.7.1 112.7.2 


wo 


OM 


Longitude of periastron at t pBT = To 


l2~7"2l 




OMDOT 


Time derivative of longitude of periastron 


rrn\ 




XOMDOT 


Non-GR contribution to derivative of longitude of periastron 


|2~731 




OM2DOT 


Second time derivative of longitude of periastron 


\'2.7.7\ 
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Table Al - continued 

Name Description Section 

Algebraic ASCII 



^4 


AO 


First binary aberration parameter 




B 


b(1 


Second binary aberration parameter 




a± 


R ACC ' 


Radial acceleration of BB-SSB 


l2.fl.2l 




PMACC 


Acceleration of BB-SSB in direction parallel to proper motion 


|2.fl.2l 


^Shk 


D_SHK 


Distance used for Shklovskii term 


12.6.21 


^AOP 


D AOP 


Distance used for annual-orbital parallax term 


|2\7.1l 




D_OP 


Distance used for orbital parallax term 


12.7.11 






Ti-th time derivative of interstellar dispersion measure 


12.6.11 


Ml 


ECC 


Eccentricity of binary orbit at t^ sr — Tq 


12.7.21 


fj 


ECCDOT 


Time derivative of eccentricity of binary orbit 


12.7.21 


i 


KI 


Binary orbital inclination angle used for Kopeikin terms 


|2. 7.112.7.21 


kt 


FDDC 


Scale of frequency-dependent delay 


|2Jjl 


M 


MTOT 


Total binary system mass 


12.7.51 


in-) 


M2 


Binary companion mass 


12.7.51 


no 


NElAU 


Mean electron density at 1 AU heliocentric radius 


I2.B.4I 


P h0 


PB 


Binary orbital period at £P sr = To 


12.7.21 




PBDOT 


Time derivative of binary orbital period 


12.7.21 


A. 


XPBDOT 


Non-GR contribution to time derivative of binary orbital period 


|2\7.5I 


T 


p, 


Shapiro delay range parameter 


I2.V.4I 


S 


S 


Shapiro delay shape parameter 


ism 


sin i 


SINI 


Sine of orbital inclination 




tu 


DMEPOCH 


Epoch of dispersion measure 


I2H 


t gn 


GLEP.n 


Epoch of n-th glitch 


EH 


£pos 


POSEPOCH 


Epoch of position 


l2.2l2.SH5.SI 




PEPOCH 


Epoch of pulse frequency 


EH 


T 


TO 


Time of periastron 


I2T2I 




TASC 


Time of ascending node, Laplace-Lagrange parameterisation 


I2T6I 




RV 


Radial velocity for Doppler correction 


I2.fl.2l2.fi.3l 


w ||Shkdot 


RVS 


Radial velocity for secular acceleration of Shklovskii term 


\2T2l 


"X 


TV 


Transverse velocity for Doppler correction 


I2TO 


a; 


Al 


Projected semi-major axis of binary orbit 


I2T2I 


i 


AlDOT 


Time derivative of projected semi-major axis of binary orbit 


I2T2I 




X2DOT 


Second time derivative of projected semi-major axis of binary orbit 


I2T7I 




SHAPMAX 


Shapiro delay parameter, z s = — ln(l — s) 


I2T21 



Tempo2: Timing model and precision 27 



Table A2. Other symbols used 
Symbol Description Section 



6c 


Unit vector in direction of increasing right ascension 


Iz.o.ol 


P 


Unit vector in direction of increasing ecliptic latitude 


12.5,31 


r 


Dopplcr shift due to BB-SSB velocity 


12.6.31 


a© 


Solar system part of timing formula 


12.51 


A A 


Atmospheric propagation delay 


12.5.21 


Aab 


Binary aberration delay 


12.7.31 


A B 


Binary part of timing formula 


12.71 


A D © 


Solar system dispersion delay 


12.5.41 


Ab© 


Solar system Einstein delay 


12.5.11 


Ab©-0 


bolar system Emstein delay lor geocentre 


12.5.11 


A 


Binary Einstein delay 


12.7.41 


A 


Einstein delay due to secular motion 


12.6.31 


A 

"FDD 


Frequency-dependent delay 


Iz.0. 11 




Atmospheric hydrostatic delay for zenith 


12.5.21 


A P 


Annual parallax delay term 


EHT51 


AlSD 


Interstellar dispersion delay 


12.6.11 


Ais 


Interstellar part of timing formula 


|2"51 


Akb 


Kopeikin delay 


I2TT1 


Arq 


Solar system Roemer delay 


|2~5"51 


A rb 


Binary Roemer delay 


I2TT1 


Ar B [| 


Radial part of binary Roemer delay 


I2TT1 


Are 


Combined binary Roemer and Einstein delay 




As© 


Solar system Shapiro delay 


|2~531 


As©2 


Second-order Solar Shapiro delay 


l2~5~5l 


A S B 


Binary Shapiro delay 


ma 


A V p 


Interstellar vacuum propagation delay 


12^21 




Atmospheric wet delay for zenith 


|2~5"21 


AL<£ N) 


Higher-order post-Newtonian correction to Einstein delay 


I^.S.ll 


Ar (A) 
jp 


Asteroid correction to Einstein delay 


Iz.o. 11 


o 


Unit vector in direction of increasing declination 


12.5.31 


V 


TT irct" T o r\ 1 o c o I c\ <tt" o ti [To rAor d m c*\ at* 
1 11 a L JjctpitlLL. IjaKl dlliiU JJcLl ellllL LLT 




e 


Source elevation angle 


ir5~2i 


6 


Orbital longitude 


I2T21 


K 


Second Laplace-Lagrange parameter 


I2T6I 


A 


Unit vector in direction of increasing ecliptic longitude 


l2~Ql 


M 


Three-dimensional proper motion vector 


ET5"51 


P 


Pulsar-Sun-observatory angle 


ES31 


<t> 


Pulse phase 


|2~S1 


<f>gn 


Contribution to pulse phase by glitch n 


E3 


f 


Geodetic latitude of the observatory 


ETS"2l 




Orbital phase-type variable for Laplace-Lagrange binary parameterisation 


l2~5~2l 




Zcroth-order approximation to $ 


I2T5I 


i>i 


Pulsar-telescope-object angle for i-th solar system body 


12331 


U! 


Longitude of periastron 


I2T21 




Angular velocity of Earth rotation 


iroi 
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Table A2 - continued 



Symbol 


Description 


Section 


A e (u) 


True anomaly of orbit as a function of eccentric anomaly 




a 


Acceleration vector of BB - SSB 


12.312.6.21 


a 


Semi-major axis of pulsar orbit 


I2T2I 


b 


Vector from BB to pulsar 


IOI2T2I 


c 


Vacuum speed of light = 2.99792458 X 10 8 m s" 1 




C 


Convenience variable in binary motion 


I2T2I 


dp 


Parallax distance = 1 AU/n 


12331 


e,eg, e r 


Orbital eccentricities 


I2T2I 


e"i , 


Basis vectors in plane of sky (i.e. a, 8 or X, /3) 


I2T2I 


D 


Interstellar dispersion constant 


123m 


DM 


Interstellar dispersion measure 


123m 


E 


Matrix relating equatorial and ecliptic coordinates 


12331 


/ 


Frequency of electromagnetic wave 


12331123m 


/S SB 


Frequency of electromagnetic wave in SSB frame 


12331123m 


G 


Gravitational constant 




H 


Height of observatory above the geoid 


12331 


i 


Orbital inclination angle 


I2T21 


k 


Displacement of BB since t SSB = i pos 


1231 12331 12331 


k 


Derivative of longitude of periastron w.r.t. true anomaly 


12331 


ko 


Constant of proportionality relating D and DM 


123m 


m% 


Binary companion mass 


12331 


ra h {Q) 


Mapping function for hydrostatic propagation delay 


12331 


m w (Q) 


Mapping function for wet propagation delay 


12331 


n 


Mean angular frequency of binary since To 


12331 


P 


Surface atmospheric pressure 


12331 


Ph 


Binary orbital period 


12331 


Q 


Frame bias, precession and nutation matrix 


12331 


R 


Earth rotation matrix 


12331 


r 


BCRS position vector of observatory 


123112331 


r @ 


BCRS position vector of geocentre 


123m 12331 


R 


BCRS displacement vector of pulsar minus observatory 


1231 


Ho 


BCRS position vector of binary barycentre at t SSB = t pos 


1231 


■Rbb 


BCRS position vector of BB at time of observation 


12331 


■Rbb 


BCRS direction unit vector of BB at time of observation 


12331 


s 


GCRS position vector of observatory 


123m 12331 


S 


Convenience variable in binary motion 


12331 


*'IS 


ISM part of geometric path length (minus value at t SSB = t pos ) 




SB 


Binary orbital part of geometric path length 


123m 


to 


Common epoch of TT, TCG and TCB 


123m 


t obs 


Terrestrial Time (~ proper time) at observatory 


1231 


t SSB 


Barycentric Coordinate Time 


1231 


t BB 


Coordinate time in frame of binary barycentre 


1231 


t psi 


Proper time in pulsar frame 


1231 


,psr 


Proper time of emission of pulse 


1231 


^framc 


Time of arrival of pulse in a given frame 


1231 


Tq 


GMq 1 c 3 , half the light travel time across the Solar Schwarzschild radius 


12331 


u m 


Gravitational potential at the geocentre due to SS bodies 


123m 


u 


Eccentric anomaly of orbit 


12331 


v® 


Velocity of the geocentre relative to the SSB 


123m 


W 


Polar motion matrix 


12331 


Wo 


Gravitational plus spin potential of Earth at the geoid 


123m 
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